function savings_ev_cubic_weighted = generate_results_cubic_external_validity(bsimax,sample,sample_wave,target_wave,tcost,covariate,winter_prevuse,SMC,savedir,fbase)
% This function calculates cost savings or energy reduction from applying
% EWM cubic rules across waves. 

% Inputs: 
% (1) bsimax: number of bootstrap runs
% (2) sample: cross-sectional dataset merged using elec_wave367_cross.csv
% and propensity scores and renamed as data_estimation_pooled_ps
% (3) sample_wave: indicates the wave the original EWM rule was estimated
% on (=3, 6, or 7)
% (4) target_wave: indicates the wave the EWM rule was applied to (=3, 6,
% or 7)
% (5) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (6) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (7) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (8) SMC: =1 use social marginal cost; =0 use retail electricity price
% (9) savedir: directory to save the savings table 
% (10) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Output: 
% (1) savings_cubic_weighted: a table summarizing the savings point estimate,
% share of households to treat, total savings for the entire state, and 95%
% confidence intervals for the point estimate 

%% Call CPLEX
addpath 'C:\Progra~1\IBM\ILOG\CPLEX_Studio129\cplex\matlab\x64_win64'

% Setup CPLEX optimization options
opt = cplexoptimset('cplex');
opt.parallel = 1;
opt.threads = 4;
opt.simplex.tolerances.feasibility = 1e-08;
opt.mip.tolerances.integrality = 1e-08;
opt.mip.strategy.variableselect = 3;
opt.mip.strategy.nodeselect = 2;
opt.mip.strategy.lbheur = 1;
opt.mip.limits.cutpasses = -1;
opt.display = 'off';

%% Gather inputs
wave=0;
[X,Y,D,n,ps,Yscale,x1_wave,prevuse_wave,Xscale] = generate_inputs_cubic_rule_reweight(sample,sample,wave,tcost,covariate,SMC);

wave=3;
sampe_wave=sample(sample.opower_paper_wave==wave,:);
[X3,Y3,D3,n3,ps3,Yscale3,x1_wave3,prevuse_wave3,Xscale3] = generate_inputs_cubic_rule_reweight(sampe_wave,sample,wave,tcost,covariate,SMC);
scaled_att_wave3= (sum(D3.*Y3.*Yscale3-((1-D3).*Y3.*Yscale3.*ps3)./(1-ps3))/sum(D3))*sum(D3)/n3;
wave=6;
sampe_wave=sample(sample.opower_paper_wave==wave,:);
[X6,Y6,D6,n6,ps6,Yscale6,x1_wave6,prevuse_wave6,Xscale6] = generate_inputs_cubic_rule_reweight(sampe_wave,sample,wave,tcost,covariate,SMC);
scaled_att_wave6 = (sum(D6.*Y6.*Yscale6-((1-D6).*Y6.*Yscale6.*ps6)./(1-ps6))/sum(D6))*sum(D6)/n6;
wave=7;
sampe_wave=sample(sample.opower_paper_wave==wave,:);
[X7,Y7,D7,n7,ps7,Yscale7,x1_wave7,prevuse_wave7,Xscale7] = generate_inputs_cubic_rule_reweight(sampe_wave,sample,wave,tcost,covariate,SMC);
scaled_att_wave7 = (sum(D7.*Y7.*Yscale7-((1-D7).*Y7.*Yscale7.*ps7)./(1-ps7))/sum(D7))*sum(D7)/n7;

X_bywave=[X3;X6;X7]; % stack three waves vertically, so that the density ratio function has data from all waves
opower_paper_wave=sample.opower_paper_wave;
Y_bywave=[Y3;Y6;Y7];
D_bywave=[D3;D6;D7];
%% Estimate g using IPW
g = Y.*(D./ps - (1-D)./(1-ps)); % elements of W(G) - W(G_0)

g3 = Y3.*(D3./ps3 - (1-D3)./(1-ps3)); 
g6 = Y6.*(D6./ps6 - (1-D6)./(1-ps6)); 
g7 = Y7.*(D7./ps7 - (1-D7)./(1-ps7)); 

g_bywave=[g3; g6; g7]; % stack g as well

%% Define X_target_wave and g_target_wave
% X and g represent the pooled sample, so we define X_target_wave, 
% g_target_wave, Yscale_target_wave for wave-specific calculation
if target_wave==3
    X_target_wave=X3;
    g_target_wave=g3;
    Yscale_target_wave=Yscale3;
    n_target_wave=n3;
    x1_target_wave=x1_wave3;
    prevuse_target_wave=prevuse_wave3;
    Xscale_target_wave=Xscale3;
    scaled_att_target_wave=scaled_att_wave3;
elseif target_wave==6
    X_target_wave=X6;
    g_target_wave=g6;
    Yscale_target_wave=Yscale6;
    n_target_wave=n6;
    x1_target_wave=x1_wave6;
    prevuse_target_wave=prevuse_wave6;
    Xscale_target_wave=Xscale6;
    scaled_att_target_wave=scaled_att_wave6;
elseif target_wave==7
    X_target_wave=X7;
    g_target_wave=g7;
    Yscale_target_wave=Yscale7;
    n_target_wave=n7;
    x1_target_wave=x1_wave7;
    prevuse_target_wave=prevuse_wave7;
    Xscale_target_wave=Xscale7;
    scaled_att_target_wave=scaled_att_wave7;
end

if sample_wave==3
    X_sample_wave=X3;
    g_sample_wave=g3;
    Yscale_sample_wave=Yscale3;
    n_sample_wave=n3;
    x1_sample_wave=x1_wave3;
    prevuse_sample_wave=prevuse_wave3;
    Xscale_sample_wave=Xscale3;
    scaled_att_sample_wave=scaled_att_wave3;
elseif sample_wave==6
    X_sample_wave=X6;
    g_sample_wave=g6;
    Yscale_sample_wave=Yscale6;
    n_sample_wave=n6;
    x1_sample_wave=x1_wave6;
    prevuse_sample_wave=prevuse_wave6;
    Xscale_sample_wave=Xscale6;
    scaled_att_sample_wave=scaled_att_wave6;
elseif sample_wave==7
    X_sample_wave=X7;
    g_sample_wave=g7;
    Yscale_sample_wave=Yscale7;
    n_sample_wave=n7;
    x1_sample_wave=x1_wave7;
    prevuse_sample_wave=prevuse_wave7;
    Xscale_sample_wave=Xscale7;
    scaled_att_sample_wave=scaled_att_wave7;
end
%% Solve for rule using cplex 
% first BS permutation- original sample
bsperm = [1:n]';
gn = g(bsperm,:);
Xn = X(bsperm,:);

% number of variables (including the constant)
k = size(X,2);
bsi=0;

tic1= tic;
% Solve for optimal treatment rule
    % generate input for cplex package
    [f,f_n,Aineq_d,Aineq_i,bineq,lb,ub,ctype,dr_hh,dr,Xu_sample_wave,gu_p,gr3,gr6,gr7,Xr3,Xr6,Xr7,Yr3,Yr6,Yr7,Dr3,Dr6,Dr7] = ...
        generate_cplex_inputs_weight(bsperm,g,X,bsi,k,gn,Xn,g_bywave,X_bywave,Y_bywave,D_bywave,opower_paper_wave,wave,sample_wave,target_wave);
    % Cost minimization with treatment decreasing in pre treatment usage
    [sol_pd, v_pd] = cplexmilp(f,Aineq_d,bineq,[],[],[],[],[],lb,ub,ctype,[],opt); 
    % with treatment increasing in pre treatment usage
    [sol_pi, v_pi] = cplexmilp(f,Aineq_i,bineq,[],[],[],[],[],lb,ub,ctype,[],opt);
    
    % get beta and v 
    if (v_pd < v_pi)
            beta = sol_pd(1:k,:);
            v    = v_pd;
        else
            beta = sol_pi(1:k,:);
            v    = v_pi;
    end
    % get solved assignment in_Ghat and cost_saving
    in_Ghat=(X_target_wave*beta>0);
    this_cost_saving=nanmean(g_target_wave.*in_Ghat)*Yscale_target_wave;

fprintf('cplex solution takes %.2f sec\n',toc(tic1))

%% Save gr, Xr for the sample wave, and calculated the scaled att for the re-weighted sample wave
if sample_wave==3
    gr_sample_wave=gr3;
    Xr_sample_wave=Xr3;
    scaled_att_sample_wave_weighted=(sum(Dr3.*Yr3.*Yscale3.*dr_hh-((1-Dr3).*Yr3.*Yscale3.*dr_hh.*ps3)./(1-ps3))/sum(Dr3))*sum(Dr3)/n3; % Remember dr_hh is sorted by X, so need to use Yr3 and Dr3 here not Y3 and D3
elseif sample_wave==6 
    gr_sample_wave=gr6;
    Xr_sample_wave=Xr6;
    scaled_att_sample_wave_weighted=(sum(Dr6.*Yr6.*Yscale6.*dr_hh-((1-Dr6).*Yr6.*Yscale6.*dr_hh.*ps6)./(1-ps6))/sum(Dr6))*sum(Dr6)/n6;
elseif sample_wave==7
    gr_sample_wave=gr7;
    Xr_sample_wave=Xr7;
    scaled_att_sample_wave_weighted=(sum(Dr7.*Yr7.*Yscale7.*dr_hh-((1-Dr7).*Yr7.*Yscale7.*dr_hh.*ps7)./(1-ps7))/sum(Dr7))*sum(Dr7)/n7;
end
%% Two-sided 95% confidence interval not estimated, set at zero for placeholder

v_pd_itr=0;
v_pi_itr=0;
v_nd_itr=0;
v_ni_itr=0;

vhats_p= -min(v_pd_itr,v_pi_itr);
vhats_n= -min(v_nd_itr,v_ni_itr);
vhats_abs = max([vhats_p vhats_n],[],2);
 
%% Export rule parameters for graphs

switch winter_prevuse % Use month_since_opower=[-21,-1] as pre-treatment months 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 % cost_savings
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end

s_wave=sprintf('%d%d',sample_wave,target_wave);

filename_coefs=sprintf('coef_ev_cubic_%s_%s%s_%s_weighted.mat',covariate,s_cost,s_winter,s_wave);
filename = sprintf('%s_%s',fbase,filename_coefs);

save(fullfile(savedir,filename),'beta',...
    'in_Ghat','x1_target_wave','prevuse_target_wave','g_target_wave','covariate','tcost','Yscale_target_wave','n_target_wave','X_target_wave','Xscale_target_wave',...
    'v','vhats_n','vhats_p','dr','dr_hh','Xu_sample_wave','X_sample_wave','gu_p','g_sample_wave','Yscale_sample_wave','Xscale_sample_wave','n_sample_wave','gr_sample_wave',...
    'Xr_sample_wave','x1_sample_wave','prevuse_sample_wave','scaled_att_target_wave','scaled_att_sample_wave','scaled_att_sample_wave_weighted');

%% Output for table 
percent=mean(in_Ghat);
savings=nanmean(g_target_wave.*in_Ghat)*Yscale_target_wave;
if (tcost)
    total_savings=savings*n_target_wave*12;
else
    total_savings=savings*n_target_wave*12/1000;
end
ci_lb=(v-prctile(vhats_abs,95))*Yscale/n_target_wave;
ci_ub=(v+prctile(vhats_abs,95))*Yscale/n_target_wave;

savings_ev_cubic_weighted=nan(1,6);
savings_ev_cubic_weighted=array2table(savings_ev_cubic_weighted,'VariableNames',{'covariate','percent','savings','totalsavings','cilb','ciub'}); 

savings_ev_cubic_weighted.covariate = {covariate};
savings_ev_cubic_weighted.percent=percent;
savings_ev_cubic_weighted.savings=savings;
savings_ev_cubic_weighted.totalsavings=total_savings;
savings_ev_cubic_weighted.cilb=ci_lb;
savings_ev_cubic_weighted.ciub=ci_ub;
end






